<<"PlotLegends`"

Simple Case

U=-(P1 δ1+P2 δ2+b12 Cos[δ1-δ2]+b13 Cos[δ1]+b23 Cos[δ2])

Plot3D[U/.{P1->0,P2->0,b12->1,b13->1,b23->1},{δ1,-π, π},{δ2,-π,π},PlotPoints->30,Mesh->True,Boxed->False,ColorFunction->(GrayLevel[((#+.1)/1.1)^(1/2)]&),AxesLabel->{"(Subscript[θ, 1])","(Subscript[θ, 2])","U"},ImageSize->Medium,PlotRange->Automatic]

Export["TorusEx1a.eps",%40,"EPS"]

Export["TorusEx1a.eps",%37,"EPS"]

Plot3D[U/.{P1->0,P2->0,b12->1,b13->1,b23->1},{δ1,-π, π},{δ2,-π,π},PlotPoints->30,Mesh->True,Boxed->False,ColorFunction->(GrayLevel[((#+3)/4)^(3)]&),AxesLabel->{"(Subscript[θ, 1])","(Subscript[θ, 2])","U"},ImageSize->Medium,PlotRange->Automatic]

Export["TorusEx1a.eps",%43,"EPS"]

ContourPlot[U/.{P1->0,P2->0,b12->1,b13->1,b23->1},{δ1,-π, π},{δ2,-π,π},PlotPoints->50,Contours->30,ColorFunction->(GrayLevel[((#+.2)/1.1)^(1/2)]&),FrameLabel->{"(Subscript[θ, 1])","(Subscript[θ, 2])"},ImageSize->Medium]

Export["TorusEx1b.eps",%60,"EPS"]

Export["TorusEx1b",%13,"EPS"]

Plot3D[U/.{P1->π/15,P2->0,b12->1,b13->1,b23->1},{δ1,-π, π},{δ2,-π,π},PlotPoints->30,Mesh->True,Boxed->False,ColorFunction->(GrayLevel[((#+3)/4)^(3)]&),AxesLabel->{"(Subscript[θ, 1])","(Subscript[θ, 2])","U"}]

Export["TorusEx2a.eps",%45,"EPS"]

ContourPlot[U/.{P1->.25,P2->0,b12->1,b13->1,b23->1},{δ1,-π, π},{δ2,-π,π},PlotPoints->50,Contours->40,ColorFunction->(GrayLevel[((#+.15)/1.1)^(1/4)]&),FrameLabel->{"(Subscript[θ, 1])","(Subscript[θ, 2])"}]

Export["TorusEx2b.eps",%58,"EPS"]

Export["TorusEx2b.eps",%47,"EPS"]

Plot3D[U/.{P1->.6,P2->-1,b12->1,b13->.5,b23->1},{δ1,-π, π},{δ2,-π,π},PlotPoints->30,Mesh->True,Boxed->False,ColorFunction->(GrayLevel[((#+3)/4)^(3)]&),AxesLabel->{"(Subscript[θ, 1])","(Subscript[θ, 2])","U"}]

ContourPlot[U/.{P1->.6,P2->-1,b12->1,b13->.5,b23->1},{δ1,-π, π},{δ2,-π,π},PlotPoints->50,Contours->50,ColorFunction->(GrayLevel[((#+3)/4)^(3)]&),FrameLabel->{"(Subscript[θ, 1])","(Subscript[θ, 2])"}]

gr1=Plot3D[U/.{P1->π/5,P2->-1,b12->1,b13->.5,b23->1},{δ1,-π, π},{δ2,-π,π},PlotPoints->30,Mesh->True,Boxed->False,ColorFunction->(GrayLevel[((#+3)/4)^(3)]&),AxesLabel->{"(Subscript[θ, 1])","(Subscript[θ, 2])","U"}]

Export["TorusEx3a.eps",%49,"EPS"]

ContourPlot[U/.{P1->π/5,P2->-1,b12->1,b13->.5,b23->1},{δ1,-π, π},{δ2,-π,π},PlotPoints->50,Contours->50,ColorFunction->(GrayLevel[((#+3.6)/4)^(3)]&),FrameLabel->{"(Subscript[θ, 1])","(Subscript[θ, 2])"}]

Export["TorusEx3b.eps",%61,"EPS"]

Export["TorusEx3b.eps",%51,"EPS"]

Equations of Motion

The system is a special case (lossless) of the 4-bus system shown in Fig. 9, Kwatny, Fischl and Nwankpa.

U=-(P1 δ1+P2 δ2+b12 Cos[δ1-δ2]+b13 Cos[δ1]+b23 Cos[δ2])
eq1=D[U,δ1]/.{P1->0,P2->0,b12->1,b13->1,b23->1}
eq2=D[U,δ2]/.{P1->0,P2->0,b12->1,b13->1,b23->1}

Reduce[eq1==0&&eq2==0&&-π≤δ1<π&&-π≤δ2<π,{δ1,δ2}]

eq1=D[U,δ1]/.{P1->0.25,P2->0,b12->1,b13->1,b23->1}
eq2=D[U,δ2]/.{P1->0.25,P2->0,b12->1,b13->1,b23->1}
Reduce[eq1==0&&eq2==0&&-π≤δ1<π&&-π≤δ2<π,{δ1,δ2}]

eq1=D[U,δ1]/.{P1->π/5,P2->-1,b12->1,b13->0.5,b23->1}
eq2=D[U,δ2]/.{P1->π/5,P2->-1,b12->1,b13->0.5,b23->1}
Reduce[eq1==0&&eq2==0&&-π≤δ1<π&&-π≤δ2<π,{δ1,δ2}]

FindRoot[{eq1==0,eq2==0},{{δ1,0},{δ2,0}}]

FindRoot[{eq1==0,eq2==0},{{δ1,2},{δ2,-2}}]

<<ProPac`

?ProPac`Dynamics`*

U

PE[δ1_,δ2_]:=-P1 δ1-P2 δ2-b13 Cos[δ1]-b12 Cos[δ1-δ2]-b23 Cos[δ2];
KE[δ1_,δ2_,ω1_,ω2_]:=ω1^2+ω2^2;
q={δ1,δ2};
p={ω1,ω2};
MakeLagrangeEquations[KE[δ1,δ2,ω1,ω2],PE[δ1,δ2],0,{0,0},t,q,p]

vars=ToExpression[Map[ToString[#]<>"[t]"&,q]]
vels=D[vars,t]

QRules=Inner[Rule,q,vars,List]
PRules=Inner[Rule,p,vels,List]

L=((T-U)/.QRules)/.PRules

Mod

U=-(P1 δ1+P2 δ2+b12 Cos[δ1-δ2]+b13 Cos[δ1]+b23 Cos[δ2])

Plot[U/.{P1->π/15,P2->0,b12->1,b13->.5,b23->1,δ2->0},{δ1,-π, 5 π},PlotPoints->500]

V=Table[U/.{P1->π/15,P2->0,b12->1,b13->.5,b23->1,δ2->0,δ1->(δ1+i 2 π)},{i,0,5}]

Plot[V,{δ1,-π, π},PlotPoints->500]

VR=Table[Rationalize[N[U/.{P1->π/15,P2->0,b12->1,b13->.5,b23->1,δ2->0,δ1->(δ1+i 2 π)}],0],{i,0,5}]

Plot[VR,{δ1,-π, π},PlotPoints->500]

TrigExpand[Series[Sin[δ1-δ2],{δ1,a,2},{δ2,b,2}]]